1 Setup

The first two sections are for reproducibility purposes. The detail of the result is shown in section 3 (Quality Control).

  • Libraries
suppressPackageStartupMessages({
    library(data.table)
    library(DESeq2)
    library(gplots)
    library(here)
    library(hyperSpec)
    library(RColorBrewer)
    library(tidyverse)
    library(VennDiagram)
    library(emoji)
    library(ggVennDiagram)
})
  • Helper files
suppressMessages({
    source(here("UPSCb-common/Rtoolbox/src/plotEnrichedTreemap.R"))
    source(here("UPSCb-common/src/R/featureSelection.R"))
    source(here("UPSCb-common/src/R/volcanoPlot.R"))
        source(here("UPSCb-common/src/R/gopher.R"))
        source(here("UPSCb-common/src/R/topGoUtilities.R"))
})
  • Graphics
pal=brewer.pal(8,"Dark2")
hpal <- colorRampPalette(c("blue","white","red"))(100)
mar <- par("mar")
  • Functions
  1. plot specific gene expression
"line_plot" <- function(dds=dds,vst=vst,gene_id=gene_id){
    #message(paste("Plotting",gene_id))
    sel <- grepl(gene_id,rownames(vst))
    stopifnot(sum(sel)==1)

    p <- ggplot(bind_cols(as.data.frame(colData(dds)),
                          data.frame(value=vst[sel,])),
                aes(x=Genotype,y=value,fill=Genotype)) +
        geom_boxplot() +
        geom_jitter(color="black") +
        scale_y_continuous(name="VST expression") + 
        ggtitle(label=paste("Expression for: ",gene_id))
    
    suppressMessages(suppressWarnings(plot(p)))
    return(NULL)
}
  1. extract the DE results. Default cutoffs are from Schurch et al., RNA, 2016
"extract_results" <- function(dds,vst,contrast,
                              padj=0.01,lfc=0.5,
                              plot=TRUE,verbose=TRUE,
                              export=TRUE,default_dir=here("data/analysis/DE"),
                              default_prefix="DE-",
                              labels=colnames(dds),
                              sample_sel=1:ncol(dds),
                              expression_cutoff=0,
                              debug=FALSE,filter=c("median",NULL),...){
    
    # get the filter
    if(!is.null(match.arg(filter))){
        filter <- rowMedians(counts(dds,normalized=TRUE))
        message("Using the median normalized counts as default, set filter=NULL to revert to using the mean")
    }
    
    # validation
    if(length(contrast)==1){
        res <- results(dds,name=contrast,filter = filter)
    } else {
        res <- results(dds,contrast=contrast,filter = filter)
    }
    
    stopifnot(length(sample_sel)==ncol(vst))
    
    if(plot){
        par(mar=c(5,5,5,5))
        volcanoPlot(res)
        par(mar=mar)
    }
    
    # a look at independent filtering
    if(plot){
        plot(metadata(res)$filterNumRej,
             type="b", ylab="number of rejections",
             xlab="quantiles of filter")
        lines(metadata(res)$lo.fit, col="red")
        abline(v=metadata(res)$filterTheta)
    }
    
    if(verbose){
        message(sprintf("The independent filtering cutoff is %s, removing %s of the data",
                        round(metadata(res)$filterThreshold,digits=5),
                        names(metadata(res)$filterThreshold)))
        
        max.theta <- metadata(res)$filterNumRej[which.max(metadata(res)$filterNumRej$numRej),"theta"]
        message(sprintf("The independent filtering maximises for %s %% of the data, corresponding to a base mean expression of %s (library-size normalised read)",
                        round(max.theta*100,digits=5),
                        round(quantile(counts(dds,normalized=TRUE),probs=max.theta),digits=5)))
    }
    
    if(plot){
        qtl.exp=quantile(counts(dds,normalized=TRUE),probs=metadata(res)$filterNumRej$theta)
        dat <- data.frame(thetas=metadata(res)$filterNumRej$theta,
                          qtl.exp=qtl.exp,
                          number.degs=sapply(lapply(qtl.exp,function(qe){
                              res$padj <= padj & abs(res$log2FoldChange) >= lfc & 
                                  ! is.na(res$padj) & res$baseMean >= qe
                          }),sum))
        if(debug){
            plot(ggplot(dat,aes(x=thetas,y=qtl.exp)) + 
                     geom_line() + geom_point() +
                     scale_x_continuous("quantiles of expression") + 
                     scale_y_continuous("base mean expression") +
                     geom_hline(yintercept=expression_cutoff,
                                linetype="dotted",col="red"))
        
            p <- ggplot(dat,aes(x=thetas,y=qtl.exp)) + 
                geom_line() + geom_point() +
                scale_x_continuous("quantiles of expression") + 
                scale_y_log10("base mean expression") + 
                geom_hline(yintercept=expression_cutoff,
                           linetype="dotted",col="red")
            suppressMessages(suppressWarnings(plot(p)))
            
            plot(ggplot(dat,aes(x=thetas,y=number.degs)) + 
                     geom_line() + geom_point() +
                     geom_hline(yintercept=dat$number.degs[1],linetype="dashed") +
                     scale_x_continuous("quantiles of expression") + 
                     scale_y_continuous("Number of DE genes"))
            
            plot(ggplot(dat,aes(x=thetas,y=number.degs[1] - number.degs),aes()) + 
                     geom_line() + geom_point() +
                     scale_x_continuous("quantiles of expression") + 
                     scale_y_continuous("Cumulative number of DE genes"))
            
            plot(ggplot(data.frame(x=dat$thetas[-1],
                                   y=diff(dat$number.degs[1] - dat$number.degs)),aes(x,y)) + 
                     geom_line() + geom_point() +
                     scale_x_continuous("quantiles of expression") + 
                     scale_y_continuous("Number of DE genes per interval"))
            
            plot(ggplot(data.frame(x=dat$qtl.exp[-1],
                                   y=diff(dat$number.degs[1] - dat$number.degs)),aes(x,y)) + 
                     geom_line() + geom_point() +
                     scale_x_continuous("base mean of expression") + 
                     scale_y_continuous("Number of DE genes per interval"))
            
            p <- ggplot(data.frame(x=dat$qtl.exp[-1],
                                   y=diff(dat$number.degs[1] - dat$number.degs)),aes(x,y)) + 
                geom_line() + geom_point() +
                scale_x_log10("base mean of expression") + 
                scale_y_continuous("Number of DE genes per interval") + 
                geom_vline(xintercept=expression_cutoff,
                           linetype="dotted",col="red")
            suppressMessages(suppressWarnings(plot(p)))
        }
    }
    
    sel <- res$padj <= padj & abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & 
        res$baseMean >= expression_cutoff
    
    if(verbose){
      message(sprintf(paste(
        ifelse(sum(sel)==1,
               "There is %s gene that is DE",
               "There are %s genes that are DE"),
        "with the following parameters: FDR <= %s, |log2FC| >= %s, base mean expression > %s"),
        sum(sel),padj,
        lfc,expression_cutoff))
    }
    
    # proceed only if there are DE genes
    if(sum(sel) > 0){
        val <- rowSums(vst[sel,sample_sel,drop=FALSE])==0
        if (sum(val) >0){
          warning(sprintf(paste(
            ifelse(sum(val)==1,
                   "There is %s DE gene that has",
                   "There are %s DE genes that have"),
            "no vst expression in the selected samples"),sum(val)))
          sel[sel][val] <- FALSE
        } 

        if(export){
            if(!dir.exists(default_dir)){
                dir.create(default_dir,showWarnings=FALSE,recursive=TRUE,mode="0771")
            }
            write.csv(res,file=file.path(default_dir,paste0(default_prefix,"results.csv")))
            write.csv(res[sel,],file.path(default_dir,paste0(default_prefix,"genes.csv")))
        }
        if(plot & sum(sel)>1){
            heatmap.2(t(scale(t(vst[sel,sample_sel]))),
                      distfun = pearson.dist,
                      hclustfun = function(X){hclust(X,method="ward.D2")},
                      trace="none",col=hpal,labRow = FALSE,
                      labCol=labels[sample_sel],...
            )
        }
    }
    return(list(all=rownames(res[sel,]),
                up=rownames(res[sel & res$log2FoldChange > 0,]),
                dn=rownames(res[sel & res$log2FoldChange < 0,])))
}
  1. extract and plot the enrichment results
extractEnrichmentResults <- function(enrichment,
                                     diff.exp=c("all","up","dn"),
                                     go.namespace=c("BP","CC","MF"),
                                     genes=NULL,export=FALSE,plot=TRUE,
                                     default_dir=here("analysis/DE"),
                                     default_prefix="DE"){
    # process args
    diff.exp <- match.arg(diff.exp)
    de <- ifelse(diff.exp=="all","none",
                 ifelse(diff.exp=="dn","down",diff.exp))

    # sanity
    if(is.null(unlist(enrichment)) | length(unlist(enrichment)) == 0){
        message("No GO enrichment for",names(enrichment))
    } else {
        # write out
        if(export){
            write_tsv(bind_rows(enrichment),
                      file=here(default_dir,
                                paste0(default_prefix,"-genes_GO-enrichment.tsv")))
            if(!is.null(genes)){
                write_tsv(
                    enrichedTermToGenes(genes=genes,terms=enrichment$id,url=url,mc.cores=16L),
                    file=here(default_dir,
                              paste0(default_prefix,"-enriched-term-to-genes.tsv"))
                )
            }
        }
        if(plot){
          gocatname <- c(BP="Biological Process",
            CC="Cellular Component",
            MF="Molecular Function")
          degname <- c(all="all DEGs",
            up="up-regulated DEGs",
            dn="down-regulated DEGs")
          lapply(names(enrichment),function(n){
            lapply(names(enrichment[[n]]),function(de){
              lapply(names(enrichment[[n]][[de]]),function(gocat){
                dat <- enrichment[[n]][[de]][[gocat]]
                if(is.null(dat)){
                  message("No GO enrichment for ",degname[de]," in category ",gocatname[gocat])
                  } else {
                    dat$GeneRatio <- dat$Significant/dat$Annotated
                    dat$Pvalue <- as.numeric(dat$FDR)
                    dat$Count <- as.numeric(dat$Significant)
                    dat <- dat[order(dat$GeneRatio),]
                    dat$Term <- factor(dat$Term, levels = unique(dat$Term))
                    ggplot(dat, aes(x =Term, y = GeneRatio, color = Pvalue, size = Count)) + 
                      geom_point() +
                      scale_color_gradient(low = "red", high = "blue") +
                      theme_bw() + 
                      ylab("GeneRatio") + 
                      xlab("") + 
                      ggtitle(paste0("GO enrichment: ",degname[de]," ",gocatname[gocat])) +
                      coord_flip()
                    }
                })
              })
            })
        }
    }
}
  • Data
load(here("analysis/salmon/dds_mergeTechRep.rda"))
dds$Genotype <- relevel(dds$Genotype,ref = "T89")

Normalisation for visualisation

vsd <- varianceStabilizingTransformation(dds,blind=FALSE)
## using 'avgTxLength' from assays(dds), correcting for library size
vst <- assay(vsd)
vst <- vst - min(vst)
#dir.create(here("analysis/DE"),showWarnings=FALSE)
#save(vst,file=here("analysis/DE/vst-aware.rda"))
#write_delim(as.data.frame(vst) %>% rownames_to_column("ID"),
#            here("analysis/DE/vst-aware.tsv"))

2 Gene of interests

2.1 Knocked-out gene

The targeted locus for line 26 is Potri.006G174000/Potra001877g14982 (Potra2n6c13821) And for the line 03, it is Potri.018G096200/Potra001273g10998 (Potra2n6c13821)

  • This is the expression level of targeted gene.
kogene <- "Potra2n6c13821"
stopifnot(kogene %in% rownames(vst))
dev.null <- line_plot(dds=dds,vst=vst,gene_id = kogene)

Even though both knockout lines have lower expression of Potra2n6c13821 than parental T89, line26 has lower expression of Potra2n6c13821 than in line3.

  • And here is preliminary statistical test to check if the expression level of the gene after knocking-out is significantly different from the parent. I used two-sided Student’s t-Test on the vst expression value
exp <- bind_cols(as.data.frame(colData(dds)),data.frame(value=vst[kogene,]))
t.test(exp$value[exp$Genotype == "line26"], exp$value[exp$Genotype == "T89"])
## 
##  Welch Two Sample t-test
## 
## data:  exp$value[exp$Genotype == "line26"] and exp$value[exp$Genotype == "T89"]
## t = -8.3332, df = 4.1565, p-value = 0.0009577
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.3056495 -0.6602609
## sample estimates:
## mean of x mean of y 
##  2.728281  3.711236
t.test(exp$value[exp$Genotype == "line3"], exp$value[exp$Genotype == "T89"])
## 
##  Welch Two Sample t-test
## 
## data:  exp$value[exp$Genotype == "line3"] and exp$value[exp$Genotype == "T89"]
## t = -4.011, df = 5.7834, p-value = 0.007582
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4919329 -0.1170351
## sample estimates:
## mean of x mean of y 
##  3.406752  3.711236

P-values for the hypothesis that the mean expression level of knocked-out lines are not equal to the parent are <0.001 for line26 and <0.01 for line3.

👉 Please note that, even though the expression levels are significantly lower, there are still read counts on the gene in knockout lines.

2.2 Other SM and LSM

suppressMessages(goi <- read_delim(here("doc/goi.txt")))
knitr::kable(goi)
Name V2
SNRPG Potra2n6c13472
SNRPB/N Potra2n11c22477
SNRPD1 Potra2n5c10994
SNRPD2 Potra2n14c27408
SNRPD3 Potra2n2c6364
SNRPF Potra2n8c17320
SNRPE Potra2n6c13821
LSM1A-B Potra2n1c1466
LSM2 Potra2n4c10338
LSM3A-B Potra2n5c11140
LSM4 Potra2n2c4537
LSM6A-B Potra2n8c17320
LSM7 Potra2n1c2324
LSM8 Potra2n1c2419
genelist <- goi$V2[!(goi$V2 == "Potra2n6c13821")]
stopifnot(all(genelist %in% rownames(vst)))
dev.null <- lapply(genelist,line_plot,dds=dds,vst=vst)

2.3 Alignment on T89 haplotypes

The reads mapped on all SM and LSM are selected and used for alignment on T89 haplotypes fasta. We do not have annotations for the haplotypes but from sequence similarity to Potra2n6c13821 we can guess where is the gene.

  • T89 Haplotype 1 (h1tg000012l) hap1

  • T89 Haplotype 2 (h2tg000021l) hap2

👉 The structure of the target gene in knockout lines is still similar to T89 parent on both haplotypes.

3 Differential Expression

dds <- DESeq(dds)
## estimating size factors
## using 'avgTxLength' from assays(dds), correcting for library size
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
  • Dispersion estimation The dispersion estimation is adequate
plotDispEsts(dds)

Check the different contrasts

resultsNames(dds)
## [1] "Intercept"              "Genotype_line26_vs_T89" "Genotype_line3_vs_T89"

3.1 Results

3.1.1 Line 26

line26 <- extract_results(dds=dds,vst=vst,contrast="Genotype_line26_vs_T89",
                                                        default_dir = here("analysis/DE"),
                                                        default_prefix = "Line26_", export=FALSE,
                                                        labels=dds$BioID,
                                                        sample_sel = dds$Genotype %in% c("line26","T89"), 
                                                        cexCol=.9, plot = TRUE, verbose = TRUE)
## Using the median normalized counts as default, set filter=NULL to revert to using the mean
## Loading required package: LSD

## The independent filtering cutoff is 0.14916, removing 28.82535% of the data
## The independent filtering maximises for 28.82535 % of the data, corresponding to a base mean expression of 0 (library-size normalised read)
## There are 8765 genes that are DE with the following parameters: FDR <= 0.01, |log2FC| >= 0.5, base mean expression > 0

3.1.2 Line 03

line03 <- extract_results(dds=dds,vst=vst,contrast="Genotype_line3_vs_T89",
                                                         default_dir = here("analysis/DE"),
                                                         default_prefix = "Line03_", export=FALSE,
                                                         labels=dds$BioID,
                                                         sample_sel = dds$Genotype %in% c("line3","T89"), 
                                                         cexCol=.9, plot = TRUE, verbose = TRUE)
## Using the median normalized counts as default, set filter=NULL to revert to using the mean

## The independent filtering cutoff is 0.14916, removing 28.82535% of the data
## The independent filtering maximises for 28.82535 % of the data, corresponding to a base mean expression of 0 (library-size normalised read)
## There are 8064 genes that are DE with the following parameters: FDR <= 0.01, |log2FC| >= 0.5, base mean expression > 0

3.2 Venn Diagram

res.list <- list("Line 26"=line26,
                                 "Line 03"=line03)

3.2.1 All DE genes

ggVennDiagram(lapply(res.list,"[[","all"), label_alpha = 0, label = "count") + 
    scale_fill_gradient(low="white",high = "dodgerblue") + 
    scale_x_continuous(expand = expansion(mult = .1))

3.2.2 Up-regulated DE genes (up in knockout)

ggVennDiagram(lapply(res.list,"[[","up"), label_alpha = 0, label = "count") + 
    scale_fill_gradient(low="white",high = "dodgerblue") + 
    scale_x_continuous(expand = expansion(mult = .1))

3.2.3 Down-regulated DE genes (down in knockout)

ggVennDiagram(lapply(res.list,"[[","dn"), label_alpha = 0, label = "count") + 
    scale_fill_gradient(low="white",high = "dodgerblue") + 
    scale_x_continuous(expand = expansion(mult = .1))

3.3 Gene Ontology enrichment

Here we used topGO to calculate GO enrichment of DEGs vs all genes. Any term having P-value less than 0.05 is shown. We didn’t calculate adjusted P-values because they did not give signicant enrichment. Please interpret the output with caution.

background <- rownames(vst)[featureSelect(vst,dds$Genotype,exp=0.00001)]
goannot <- prepAnnot(mapping = "/mnt/picea/storage/reference/Populus-tremula/v2.2/gopher/gene_to_go.tsv")

3.3.1 DEGs from Line26

res.list <- list("Line 26"=line26)

suppressMessages(enr.list <- lapply(res.list,function(r){
  lapply(r,topGO,background=background,annotation=goannot,alpha=0.05,p.adjust="none")
}))

suppressWarnings(extractEnrichmentResults(enr.list))
## [[1]]
## [[1]][[1]]
## [[1]][[1]][[1]]

## 
## [[1]][[1]][[2]]

## 
## [[1]][[1]][[3]]

## 
## 
## [[1]][[2]]
## [[1]][[2]][[1]]

## 
## [[1]][[2]][[2]]

## 
## [[1]][[2]][[3]]

## 
## 
## [[1]][[3]]
## [[1]][[3]][[1]]

## 
## [[1]][[3]][[2]]

## 
## [[1]][[3]][[3]]

3.3.2 DEGs from Line03

res.list <- list("Line 03"=line03)

suppressMessages(enr.list <- lapply(res.list,function(r){
  lapply(r,topGO,background=background,annotation=goannot,alpha=0.05,p.adjust="none")
}))

suppressWarnings(extractEnrichmentResults(enr.list))
## [[1]]
## [[1]][[1]]
## [[1]][[1]][[1]]

## 
## [[1]][[1]][[2]]

## 
## [[1]][[1]][[3]]

## 
## 
## [[1]][[2]]
## [[1]][[2]][[1]]

## 
## [[1]][[2]][[2]]

## 
## [[1]][[2]][[3]]

## 
## 
## [[1]][[3]]
## [[1]][[3]][[1]]

## 
## [[1]][[3]][[2]]

## 
## [[1]][[3]][[3]]

3.3.3 Common DEGs between two lines

res.list <- list("CommonLine26andLine03"=list("all"=intersect(line26$all,line03$all),
                                                                                            "up"=intersect(line26$up,line03$up),
                                                                                            "dn"=intersect(line26$dn,line03$dn)))

suppressMessages(enr.list <- lapply(res.list,function(r){
  lapply(r,topGO,background=background,annotation=goannot,alpha=0.05,p.adjust="none")
}))

suppressWarnings(extractEnrichmentResults(enr.list))
## [[1]]
## [[1]][[1]]
## [[1]][[1]][[1]]

## 
## [[1]][[1]][[2]]

## 
## [[1]][[1]][[3]]

## 
## 
## [[1]][[2]]
## [[1]][[2]][[1]]

## 
## [[1]][[2]][[2]]

## 
## [[1]][[2]][[3]]

## 
## 
## [[1]][[3]]
## [[1]][[3]][[1]]

## 
## [[1]][[3]][[2]]

## 
## [[1]][[3]][[3]]

4 Session Info

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] LSD_4.1-0                   topGO_2.54.0               
##  [3] SparseM_1.81                GO.db_3.18.0               
##  [5] AnnotationDbi_1.64.1        graph_1.80.0               
##  [7] limma_3.58.1                ggVennDiagram_1.4.9        
##  [9] emoji_15.0                  VennDiagram_1.7.3          
## [11] futile.logger_1.4.3         lubridate_1.9.3            
## [13] forcats_1.0.0               stringr_1.5.1              
## [15] dplyr_1.1.4                 purrr_1.0.2                
## [17] readr_2.1.4                 tidyr_1.3.0                
## [19] tibble_3.2.1                tidyverse_2.0.0            
## [21] RColorBrewer_1.1-3          hyperSpec_0.100.0          
## [23] xml2_1.3.6                  ggplot2_3.4.4              
## [25] lattice_0.21-8              here_1.0.1                 
## [27] gplots_3.1.3                DESeq2_1.42.0              
## [29] SummarizedExperiment_1.32.0 Biobase_2.62.0             
## [31] MatrixGenerics_1.14.0       matrixStats_1.2.0          
## [33] GenomicRanges_1.54.1        GenomeInfoDb_1.38.5        
## [35] IRanges_2.36.0              S4Vectors_0.40.2           
## [37] BiocGenerics_0.48.1         data.table_1.14.10         
## 
## loaded via a namespace (and not attached):
##  [1] DBI_1.2.0               bitops_1.0-7            deldir_2.0-2           
##  [4] formatR_1.14            testthat_3.2.1          rlang_1.1.2            
##  [7] magrittr_2.0.3          RSQLite_2.3.4           compiler_4.3.1         
## [10] png_0.1-8               vctrs_0.6.5             pkgconfig_2.0.3        
## [13] crayon_1.5.2            fastmap_1.1.1           XVector_0.42.0         
## [16] labeling_0.4.3          caTools_1.18.2          utf8_1.2.4             
## [19] rmarkdown_2.25          tzdb_0.4.0              bit_4.0.5              
## [22] xfun_0.41               zlibbioc_1.48.0         cachem_1.0.8           
## [25] jsonlite_1.8.8          blob_1.2.4              highr_0.10             
## [28] DelayedArray_0.28.0     BiocParallel_1.36.0     jpeg_0.1-10            
## [31] parallel_4.3.1          R6_2.5.1                bslib_0.6.1            
## [34] stringi_1.8.3           brio_1.1.4              jquerylib_0.1.4        
## [37] Rcpp_1.0.11             knitr_1.45              Matrix_1.6-0           
## [40] timechange_0.2.0        tidyselect_1.2.0        rstudioapi_0.15.0      
## [43] abind_1.4-5             yaml_2.3.8              codetools_0.2-19       
## [46] KEGGREST_1.42.0         withr_2.5.2             evaluate_0.23          
## [49] lambda.r_1.2.4          Biostrings_2.70.1       pillar_1.9.0           
## [52] KernSmooth_2.23-22      generics_0.1.3          vroom_1.6.4            
## [55] rprojroot_2.0.4         RCurl_1.98-1.13         hms_1.1.3              
## [58] munsell_0.5.0           scales_1.3.0            gtools_3.9.5           
## [61] glue_1.6.2              lazyeval_0.2.2          tools_4.3.1            
## [64] interp_1.1-5            locfit_1.5-9.8          latticeExtra_0.6-30    
## [67] colorspace_2.1-0        GenomeInfoDbData_1.2.11 cli_3.6.2              
## [70] futile.options_1.0.1    fansi_1.0.6             S4Arrays_1.2.0         
## [73] gtable_0.3.4            sass_0.4.8              digest_0.6.33          
## [76] SparseArray_1.2.2       farver_2.1.1            memoise_2.0.1          
## [79] htmltools_0.5.7         lifecycle_1.0.4         httr_1.4.7             
## [82] statmod_1.5.0           bit64_4.0.5
 

drawing

Created by Fai Theerarat Kochakarn

theerarat.kochakarn@umu.se